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FOREWORD 


The research reported herein was supported by NASA Contract NAS8-29584. 
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Laboratory, Marshall Space Flight Center, was the scientific monitor, and 
support was provided by Messrs. John Enders and Harry Verstynen of the 
Aviation Safety Technology Branch, Office of Advanced Research and Tech- 
nology (OAST), NASA Headquarters. 

The research reported in this document is concerned with the results of 
a study of atmospheric boundary layer flow over surface obstacles using 
an advance turbulent closure scheme involving the turbulence kinetic energy 
equation. It is part of our investigation to determine the wind environments 
that are induced by buildings, etc., in the vicinity of airports, heliports, 
STOL-ports, etc., relative to assessing and defining the possible hazardous 
flying conditions which could result from these wind environments. The last 
chapter of the report (Chapter IX) represents an analysis of the implications 
of the calculated induced wind environments on aeronautical safety. 
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turbulence kinetic energy 
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Dissipation term of turbulence kinetic energy 
equation 

C Stretched vertical coordinate 
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X Length of separation bubble 

V Kinematic viscosity 

p Density 

a Root mean square value of fluctuating u velocity 

component 

Root mean square value of fluctuating v velocity 
component 

Root mean square value of fluctuating w velocity 
component 

T Turbulent shear stress 

u Convection velocity 

^ Vertical coordinate stretching function 
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CHAPTER I 


INTRODUCTION 

The flow fields induced around buildings: and .other 
structures by the surface winds are of considerable impor- 
tance with respect to the design of aeronautical facilities, 
particularly in the case of helicopter and V/STOL facilities 
in metropolitan areas where these aircraft may be operated 
from the tops of buildings. Severe updrafts and downdrafts 
and regions with high turbulence levels make take-off and 
landing of aircraft near buildings hazardous. These hazards 
can be attenuated by modifying airport design or aircraft 
operation procedures so as to reduce the influence of the 
flow fields caused by the wind. This, however, requires 
methods for obtaining quantitative descriptions of the 
atmospheric flow at the aerodrome. 

One method is to measure the flow conditions at 
existing aerodromes. This method is impractical because of 
the large number of measurements required and the vast 
amount of data which would have to be reduced. 

The above comments are not meant to suggest that wind 
measurements at airports are not required, but that entire 
flow field patterns cannot practically be established in 
this manner. Measurements with one or two individual masts 
could, however, be supplemented with analytical methods 
which would allow logical extrapolation of these measure- 
ments and thus allow predictions of the overall wind pattern 



about the airfield to be made. McManus [1]^ observed that 
neither time histories nor rms values measured by aircraft 
agree with values determined from an individual mast and 
concluded that the turbulence and mean wind pattern seen by 
the aircraft was defined by the geometry of the airfield. 

Since geometry plays a significant role, it is further 
obvious that rather extensive measurements made at one air- 
field could not be employed to design few facilities whereas 
numerical models would have this flexibility. Moreover, 
analytical models of winds about irregular surface geometries 
are also of interest to airplane design and flight simulation 
(see [2 , 3] ) . 

The need of an analytical approach is thus indicated. 

A solution of the complete equations of motion for turbulent 
flow would provide the desired information but such solu- 
tions are not practical because of the great expense of com- 
puter time required to achieve results. 

The concept of a disturbed boundary layer with turbu- 
lence modeled after the Prandtl mixing length hypothesis was 
used by Frost, Maus, and Simpson [4] to approximate a 
solution of the general equations. The two-dimensional 
boundary layer equations were solved for flow over semi- 
elliptical cylinders. 

The significant conclusions of the above mentioned 

^Numbers in brackets refer to similarly numbered 
references in the list of references. 
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study were ; 

1. A localized maximiam in wind speed is produced at 
the top of a semi-elliptical body, the maximum velocity 
increasing with decreasing elliptical aspect ration. 

2. Increased surface roughness decreases the wind 
velocity in the boundary layer. 

3. Decreasing elliptical aspect ratio and/or 
increasing surface roughness cause a larger separation 
region upstream of the body. 

The question of the validity of the Prandtl mixing 
length theory for atmospheric flow in disturbed regions 
gives a degree of uncertainty to the above conclusions, how- 
ever. Moreover, this method of solution does not give any 
information about the turbulence structure of the flow field. 

The purpose of this report is to develop solutions 
which will shed light on the validity of the results obtained 
from the Prandtl mixing length hypothesis. Flow fields are 
calculated with the boundary layer equations using equivalent 
boundary conditions employed in [4] , but with the equation 
of conservation of turbulence kinetic energy used to model 
the turbulence kinetic energy statistics. If the flow 
fields predicted' by the two methods of solution do not 
differ significantly, then both models can be employed with 
greater confidence. Of the two methods of solutions, the 
Prandtl mixing length method has the advantage of requiring 
less computation time, whereas, the turbulence kinetic 
energy method has the advantage of allowing for the 
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turbulence to feed from the mean flow via the energy equa- 
tion giving physically more meaningful values of turbulence 
intensity levels. 

The method developed in this study can be applied to 
flow over very general two-dimensional bodies, however, for 
the purposes of comparison with the method of Frost, Maus, 
and Simpson [4] , the solutions presented will be limited to 
semi-elliptical cylinders with various aspect ratios and 
surface roughnesses. 



CHAPTER II 


TURBULENCE KINETIC ENERGY MODELS 


The equation of conservation of turbulence kinetic 
energy is employed to close the equations of continuity and 
conservation of momentum. In order to achieve closure, a 
relationship between the turbulence kinetic energy and the 
turbulent shear stress is hypothesized. Byrne [5] , Lee and 
Harsha [7] , and Bradshaw, Ferris, and Atwell [7] have pro- 
posed the relationship. 
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9u* 

9z* 



( 1 ) 


where t is the turbulent shear stress. 


T = -pu'w' 


e* is the turbulence kinetic energy. 


e* 



V* 2 + w' 2 ) 
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and aj is a "universal" constant. 

For a neutral atmosphere, it has been experimentally 
determined, as reported by Fichtl [5] , that: 
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where 




= ,V^)V2 


and 


n* = /FTp 


Therefore, 


i* = 


|-(u' 2 + V' 2 + 


= iu*2 

o 'J- J. 


[!u]' . f!vi= ^ rfwi 

M [u*j 


=: 6 


(3) 


This result is strictly valid for a horizontally homogeneous 
neutrally stratified boundary layer; however, we shall 
assume that it also hold for turbulent flow over an obstruc- 
tion. Thus, from Equation 1, for a neutral atmosphere, 
a^ 1/6. 

The steady state turbulence kinetic energy equation 
is (in tensor notation) from Harsha [9] : 
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1 3x| 


3u* 

= -u!u'. •= — I- - V 

ID 3x| 


Convection Production Dissipation 


- '“i'E''''’ * 

Diffusion 


The production term represents the production of 
turbulence kinetic energy by the working of the flow against 

the turbulent stresses t = -p u ! u . 

ID 

The diffusion term is composed of the gain of energy 
by flow down pressure gradients, -3 (u| p/p) /3x| , the trans- 
port of turbulence energy by large eddies, -3 (uT~"e^li/3x^ , 
and a gradient-diffusion term v3^e*/3xt^. If the gradient 
diffusion term is considered to be negligible compared to 
the convective diffusion term, the resulting turbulence 
energy budget, with the continuity and momentum equations, 
gives a hyperbolic system of equations. Bradshaw, Ferris, 
and Atwell [7] neglect the gradient diffusion term and use 
as a diffusion term, 3(aje* u)/3z* where u is a convection 
velocity defined by. 


G(z*/6) 


with G modeled from experimental data as; 


G = w' ^ + w' e*' 


7 


where t is the maximum shear stress across the velocity 
profile and S is the boundary layer thickness. The resulting 
hyperbolic system of three equations is solved by the method 
of characteristics. 

While the solution of a hyperbolic system of three 

equations may require less computer time than a parabolic 

system, the complexity of the hyperbolic solution increases 

greatly if additional equations such as species or total 

energy equations are added to the system. References [9, 6, 

f k 9e*l 

5] incorporate a gradient diffusion term, 3 /3z* 

where k is an effective viscosity defined by. 



The symbol Sc^ is an effective Schmidt number which 
represents the ratio of turbulent eddy viscosity to a dif- 
fusion coefficient for energy. 

The resulting turbulence kinetic energy equation is 
parabolic and can be solved simultaneously with the momentum 
and continuity equations . 

This latter approach is also applicable in this study 
since the primary interest lies in the lower portion of the 
atmospheric boundary layer where the average eddy size is 
small compared to that in the upper region of the boundary 
layer, and hence the magnitude of the convective diffusion 
term is small compared to the production and dissipation 
terms [10] . Consequently, the gradient diffusion term is 
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used in order to simplify the solution of the equations. 

An approximation must be made for the dissipation 
term, e = v ( 9u^y3x|) , in order to obtain a form suitable 

for computation. Byrne and Lee [11] propose a dissipation 
term, 

e = 

where 6 is the boundary layer thickness, and a^ is a param- 
eter dependent upon the location of the maximum shear stress, 

T . This dissipation term, however, does not satisfy the 
max 

condition of constant shear stress which exists for atmo- 
spheric boundary layers over flat terrain. 

Bradshaw and Ferris [7] propose that, 

e = u*^L (6) 

where, 

uj = /t/p = (a^e*)^^^ 
and 


L 


< 2 *, near the surface 

0.095y^^^f in outer layer 


where k is the von Karman constant, and y is the inter- 
mittency . 

It can be shown that for steady flow with no pressure 
gradient, corresponding to a logarithmic velocity profile 
where production equals dissipation, the dissipation term is 
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given by. 


^ k(z* + z*) 

where z* is the surface roughness height. Thus, e satisfies 
the initial upstream condition of a logarithmic approaching 
velocity profile. For z* >> z*, the value of e given by (7) 
is approximately equal to that given by (6) . Therefore (7) 
has been chosen for the dissipation function in this report. 
This dissipation function is valid far upstream, but whether 
it applies in the vicinity of the body is not known and must 
eventually be established by experiment. 

The Governing Equations 

The equations of motion for steady, incompressible 
flow within the atmospheric boundary layer are well estab- 
lished [4, 12]. The continuity equation is; 


9u* , 3w* _ 
9x* 9z* 


( 8 ) 


The longitudinal momentum equation is; 


u 


* iHl + „* 


9x* 


9z* 


= -1/p 


IP- 

9x’' 


9z* 


T + fv* 
xz 


(9) 


where f is the Coriolis parameter. 

Below a height of about 50 m the apparent force pro- 
duced by the Coriolis effect due to the Earth's rotation is 
negligible compared to the vertical gradient in t [4] , so 
the term fv* can be omitted from (9) . 


10 



Based on the approximations made for the diffusion 
and dissipation terms, the turbulence kinetic energy equa- 
tion for steady, two-dimensional flow is. 


u 


* 9e* . 9e^ ^ 3 

8x* 9z* 3z* 


3e* 


Sc 3z* 
e 


+ X 


3u* _ ui^ 

p 3z* < (z* + z*) . 


( 10 ) 


Substituting Equations (1) , (3) , and (5) into (9) and 

(10) , the governing equations become: 


Continuity, 

3u* 3w* 

3x* 3z* 


( 11 ) 


Momentum, 


u* 


3u* 

3x* 


+ w* 


3u* 

3z* 


- 1 3p 3 f aie* 9u*' 
p 3x* 3z* [ I 3u*/3z* I 3z*^ 


( 12 ) 


Turbulence Kinetic Energy, 


* 3e* , * 3e* 

_ 3 r aie* 

3e*] 

aie* 

'3u*'| 

^ 3x* 3z* 

3z* [Sc^l 3u*/9z* 1 

3z*J 

1 3u*/3z* 1 

[3z*J 


(aie*) 
K (Z* + 


3/2 


(13) 
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CHAPTER III 


NUMERICAL SOLUTION OF THE GOVERNING EQUATIONS 

The continuity, momentum, and turbulence kinetic 

energy equations form a closed set of nonlinear, parabolic 

partial differential equations given by Equations (11) , 

(12), and (13) . The boundary conditions for these equations 

are: u* = 0, e* = uJ^/aj, at z* = 0, and u* = 

e* = E^, at z* = where E^ is the turbulence kinetic 

energy corresponding to a logarithmic velocity profile, 

is a height above which the influence of the body is assumed 

negligible, and u? (z* ) is the velocity at height z* 

log max max 

corresponding to a logarithmic velocity profile. The equa- 
tions are solved numerically by the procedure described 
below. 

The pressure gradient term in Equation (12) is approxi- 
mated by the pressure distribution of the inviscid solution 
for flow over an ellipse [4] . The pressure distribution 
used corresponds to the streamline nearest to the body which 
allows stable solutions of the momentum and turbulence 
kinetic energy equations to the top of the ellipse. This 
"non-separating" streamline is determined by iteratively 
introducing streamlines further and further from the body 
until a pressure distribution which permits a complete, 
stable solution is found. Jackson and Hunt [13] have, 
through a perturbation analysis, analytically confirmed the 
validity of this pressure distribution for the outer layer 
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of the flow field over a surface hump. (See also Reference 
[4].) From Bernoulli's equation: 


1 dp 
p dx* 



(14) 


the pressure gradient term in Equation (12) can be expressed 
in terms of U|, which is the velocity given by the potential 
flow solution along the non-separating st-reamline. 

Adopting a characteristic length b, equal to the 
height of the ellipse, a characteristic velocity U^, equal 
to (u^/<) In (3b + z*)/z* for the initial velocity profile, 
and a characteristic turbulence kinetic energy E^, equal to 
u^/a^ for the initial turbulence kinetic energy profile, the 
governing equations can be expressed in the following non- 
dimensional form. 


9uaw_- 


(15) 


9u , 9u 

U + W 

9x 9 z 


dU a E „ 

U — ^ + -L-£ 1- 
“^e dx U2 9z 


e 9u 
9u/9 z I 9z 


(16) 
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— + w — = 

SL E 

10 9 

^ e 

9e' 

1 -L 

9u 

9x ^ " 9z 

Sc U2 9z 
e °° 

L |lu7aTl 

’ 9zJ 

1 ^ 



^3/2-1511/2 

®o e'^2 


K U, 


(z + z^) 


(17) 


where : 


e = e*/E^ 
u = u*/U^ 


13 



w = w*/U^ 

= U*/U 
e e ' 

X = x*/b 

z = z*/b 

Zq = z*/b 


In the finite difference scheme used to solve the 
governing equations, it is advantageous to transform the 
vertical coordinate in order to compress the vertical length 
scale near the surface such that more grid points will occur 
in the region where the velocity gradients are most severe. 
The transformation is in the form of a stretching function 
C = f(z). Applying the transformation to the vertical 
derivative results in the relationship: 


9z 


f ’ (z) 


9 

9 ^ 


Defining the function. 


$(C) s f • (z) = f ' [g(0 ] 


the governing equations becomes 


Continuity, 


= 0 

9x 9C 


(18) 


Mom.entum, 
u + w$ 


9u 
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e dx 


a E 
1 o 
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9 f $ e 9u' 
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Turbulence Kinetic Energy, 
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( 20 ) 


Finite Difference Form 

The preceding equations are solved by an implicit, 
finite difference technique. The flow field over the 
cylinder is assumed to be represented by Cartesian grid 
system shown in Figure 1. Curvature effects are assumed 
negligible since the flow field within the stagnation 
regions has already been neglected by the assumption that 
the pressure distribution imposed on the flow is given by 
that along the non-separating streamline [4] . The equations 
of motion can be solved by approximating the derivatives 
with finite difference forms expressed in terms of nodal 
points of the grid. This approximation results in a set of 
equations involving the unknown values at an x-station m+1 
and known values at x-station m. The derivatives of the u 
velocity component are expressed as: 


3u 

3x 


m+l,n 


^m+l,n ^m,n 
Ax 


3u 

3? 


m+1 , n 


'^m+l,n+l '^m+l,n-l 

2A? 
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0 ) 



(b) ASSUMED NUMERICAL COORDINATES 


Figure 1. 


Coordinate System Over the Elliptical Cylinder 



9"u 


m+l,n 


^m+l,n+l ^^m+l,n ^m+l,n-l 

( A ?)2 


The derivatives of the other variables are expressed simi- 
larly, Since the momentum and turbulence kinetic energy 
equations are nonlinear, an iterative procedure must be 
introduced into the numerical scheme - 

The finite difference equations are made linear by 
replacing the dependent variables which appear to powers 
other than the first by a factored product of the variable 
to the first power and the remainder which is approximated 
with values of the variable from the previous iteration. 

For the first iteration at an x-station m+1, the known 
values of the variable (u or e) at station m are used to 
linearize the equations. The system is solved at m+1, and 
this solution is used to linearize the equations for the 
second iteration. This procedure is repeated until the 
desired degree of convergence is attained. The method is 
then repeated at m+2 using the solution at m+1 to initially 
linearize the equations. In this way the numerical pro- 
cedure marches downstream. 

The solution procedure is initiated at a prescribed 
distance upstream with a logarithmic velocity profile. The 
initial turbulence kinetic energy profile is obtained by 
iterating on e with the logarithmic u profile fixed. The 
solution procedure then advances to the first x-station. 

The first step in the iteration loop is to solve the 
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III Ml nil nil 
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momentiam equation for u, using values of e obtained from the 
previous iteration, or in the case of the first iteration at 
m+1, from the solution at m. The momentum equation. Equation 
(19), is repeated here for convenience. 


9u . . . 9u ,, 
“ as 3? ' " 


dU 
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+ — ' , <j) - — 

e dx U^_ 


^ e 9u 

TFT |Su/9c| 9C 


Note that the last term on the right side can be written. 


^1^0 1^1 , 9e 
1^1 ^ 9^ 


where , 


- _ 9u/9c 

|9u/9? 


+1, 9u/9c > 0 
-1, 9u/9c < 0 


Since 


4> ^ [+1, 4> > 0 

JWl [-1, $ < 0 


$ — 
9^ 


,T^ «) = l»l 


9e 


The finite difference form of the momentum equation is then. 
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'^m+ljn Ax ^n'^m+l,n 


dU 


U 
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2Ac 

T _ T 

^m+l,n+l ^m+1, n-1 
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(21) 


where the superscript T denotes a trial value from the pre- 
vious iteration. 
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Rearranging terms, we obtain. 
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which has the form. 


Au., .t+Bu-, +Cu., -=D, (2<n<N-l) 

n m+l,n+l n m+l,n n m+l,n-l n 


(23) 


where N is the number of nodal points at m+1, and 


A = 
n 


w 

n m+l/n Ax 

2 Ac 


B = u 


n 


C = 
n 


m+l,n 

-A 

n 


D = 
n 


U 


e dx 
T 


3. E 

Ax + A|$ I 

m+1 ” 


T _ T 

^m+l,n+l ^m+l,n-l 


Ax 

AC 


^ '^m+l,n ^^m,n 


Applying these equations at each of the nodal points 
between the surface and outer boundary results in a system 
of N - 2 algebraic equations which yield a tridiagonal 
coefficient matrix which can be solved by the following 
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recursion formula: 


u., =E u,- ,t+F 
m+l,n n m+l,n+l n 


and. 


'^m+l,n-l ^n-1 ^m+l,n ^n-1 


( 24 ) 


Substituting into Equation (23) , 


A u ,, ,, + B u ,, + C [E tU ,, + F ,] = D 

n m+l,n+l n m+l,n n n-1 m+l,n n-1 n 


and solving for 


-A u - D - C F T 

_ n m+l,n+l n n n-1 

Wl,n B+CE, B+CE, 
n n n-1 n n n-1 


which means. 


E IL_ 

n B + C E 


n 


*1 

n n-1 


D - C F T 
p _ n n n-1 

""n B^ + C E^ T 
n n n-1 


(25) 


(26) 


Boundary Conditions on u 

The inner boundary condition of no slip at the ground 

in Equation (24) gives Ej and F^ as zero. Equations (25) 

and (26) are used to calculate the E 's and F 's from n = 2 

n n 

outward to n = N - 1 . The outer boundary condition. 
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u. 


Wl,N 


K 


In 


6.3+ z 


which was chosen based on computational experience with the 

program, is used in Equation (24) to begin calculation of 

the u ' s from n = N - 1 inward to n = 2 . 
n 

The velocity profile in the vertical direction is 
solved by integrating the continuity equation. Equation (18), 


w(?) 



Approximating the integral by the trapezoidal rule, the 
vertical velocity is calculated outward from the surface 
with. 


Ag 

'^m+l,n ''^m+l,n-l 2 


1 u . T , - u . 

1 m+l,n-l m, n-1 


n-1 


Ax 


1 ^m+l,n ^m,n 


n 


Ax 


(27) 


We now direct attention to the turbulence kinetic 
energy equation. Equation (20) , which is repeated here for 
convenience. 


9 e , . 3 e 

u + $w 
8x 


^i®o . 3 f e 3e' 

|3u/3c| 


+ ai 1$ I 


e 


8u 


-3/2 Tjt 1 / 2 

^o 
K U 


, 3/2 


(z + z^) 
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The finite difference form of the equation is 


T ®m+l,n ®m,n . * ,.T ®m+l,n+l ®m+l,n-l 

"m+l,n * *n "m+l,n 24c 


a,E 
1 o 


1 ^ 


Sc U2 I^n' ! 
0 00 . 


>• ' 

®m+l,n+l ®m+l,n-l 

2 

2AC 1 

2AC J 


T 

%+l,n+l 

T 

'^m+l,n-l 


+ e 


m+l,n 


^m+l/n+1 ^^m+I,n ^ ^m+l,n-l 


(AC) 


2AC 

— ^ 

®m+l,n+l ®m+l,n-l 

T _ T 

'^m+l,n+l '^m+l,n-l 

m+1 , n 

2AC J 


2Ac 


^m+l,n+l '^m+l,n— 1 


^^m+l,n+l ^m+l,n ^ 


(AC) 


l^m+l,n ^m+l,n-l 


(AC) 

T _ T 

'^m+ljn+l '^m+l,n-l 


+ ^il'^nl ®m+l,n 


2AC 


m 1 / 2 

3 / 2 pi / 2 J- 

1 ^o ^m+l,n m+l,n 


K U 


(z + z^) 


(28) 


Rearranging terms. Equation (28) can be written as. 


G e. 


+ H e 


+ I e 


n m+l,n+l n m+l,n n m+l,n-l n 


= J, (2<n<N-l) 


(29) 


where , 
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$ a E I $ I e"^ - 20 ”^ + e"^ 

_ n m+l,n Ax _ 10 * ^ ^m+l,n+l m+l,n ^m+l,n-l 


AC SCg 


(2AC) 


2AC 

®m+l , n 

2AC 

T T 

%+l,n+l " %+l,n-l 

2 

(AC)"' 

T _ T 

''‘m+1, n+1 ~ ''‘m+l,n-l 


'm+1 , n 
2 


'^m+l/n+1 '^m+l/n— 1 


'^m+l,n+l ^m+1 , n ^ 


I T _ T I 

I ^m+lyii ~ ^m+l,n-l^ 


H = 
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m+1, n 2Ac 
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I . I m+l/n+1 ^m+l/ii— 1 . 

- 'm ** 
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a E U> 
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®m+l,n-l 
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J = e u 
n m,n m+l,n 


As with Equations (23) , Equations (29) form a system of 
N - 2 equations which yield a tridiagonal coefficient matrix 
which can be solved by the formula. 


^m+l,n-l ^n-1 ®m+l,n ^n-1 


(30) 


where , 


n 


^n H + I K 


(31) 


n 


n n-1 


L„ = 


J - I L T 
n n n-1 


n H + I K T 
n n n-1 


(32) 


Boundary Conditions on e 

The inner boundary condition on e is. 


l"wl 


'm+1,1 p E^ 


where , 


"*m+l,l 


%+l,2 

ln[ (Z 2 + z^) ] 


(33) 


which assumes that the logarithmic law is obeyed at the 
wall. Substituting into Equation (30) gives, 

Kj = 0 
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Ll = 




Equations .( 31) and (32) are used to calculate the K^^'s and 
L^'s outward from n=2ton=N-l. The outer boundary 
value of e = 1.0 is that for an undisturbed boundary layer 
at the same height above the surface. This boundary condi- 
tion is used in Equation (30) to begin calculation of the 

e ' s inward from n = N to n = 1 . 
n 

In summary, the procedure for calculating the u and w 
velocity profiles and turbulence kinetic energy profile at 
m+1 is: 

1. Calculate A , B , C , and D from known values. 

n n n n 

2. Calculate E and F by Equations (25) and (26) 

n n 

from n=lton=N-l, starting with Ej = Fj = 0. 

3. Calculate U ,, inward from the outer boundary 

m+ 1 , n 

condition to n = 2 from Equation (23) . 

4. Calculate w ,, in steps from the inner boundary 

m+ 1 , n 

condition Wi = 0 to the outer edge by Equation 
(27) . 

5. Calculate G , H , I , J from prescribed values. 

n n n n 

6. Calculate K and L by Equations (31) and (32) 

n n 

from n = lton = N- 1, starting with == 0, 

L, = (u* T U )^/a,E . 

1 *m+l,l oo' ' 1 o 

7. Calculate e inward from n = N to n = 1 by 

m+ 1 f n 

Equation (30) with ^ i^o 

undisturbed boundary layer. 
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8. Return to step 1 using newly calculated u, w, and 
e profiles, and repeat procedure until desired 
convergence is attained. 


Convergence 

The convergence of the numerical scheme given above 
is strongly affected by the magnitude of the pressure 
gradient. In order to minimize computational time, a 
variable step size in the x-direction, was incorporated. 
The step size is adjusted such that larger values of Ax 
occur in regions of small pressure gradient and smaller 
values in regions of large pressure gradient. A step size 
determined by. 


Ax « 


dU 


U 


e dx 


- 1/2 


has been found to give satisfactory results. 

Convergence criteria established in Reference [4] for 
grid spacing Ax, and AC to assure convergence of the itera- 
tion process and to minimize truncation error were also 
observed . 
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CHAPTER IV 


RESULTS AND DISCUSSION 

Niamerical solutions of the governing equations have 
been carried out for elliptical cylinders with aspect ratios 
(length of semi-major axis divided by the height of cylinder) 
of 2/1, 4/1, and 10/1. The results of these calculations 
and comparison with a previous solution of the conttinuity 
and momentum equations based on an eddy viscosity model [4] 
are presented. 

The calculation procedure described previously is 
initiated with the logarithmic velocity profile character- 
istic of the approaching wind at a distance sufficiently far 
upstream of the body that the pressure disturbance created 
by the body is negligible compared to the free stream 
pressure. Marching downstream, the velocity profiles are 
calculated using u = u*^^ln[(6.3 + z^)/z^]/k as the upper 
bounding condition. 

The parameters affecting the solution are the aspect 
ratio of the elliptical body, a, and the surface roughness 
of the upstream terrain, z^, (z^ is assumed the same for the 

body) . 
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CHAPTER V 


COMPARISON OF TURBULENCE KINETIC ENERGY MODEL 
WITH PRANDTL MIXING LENGTH MODEL 


Velocity profiles calculated by the turbulence 
kinetic energy (TKE) model with aj=l/6 and the Prandtl 
mixing length (PML) model of [4], are shown in Figures 2, 3, 
and 4. The flow velocities calculated by the TKE model lag 
those calculated by the PML model in the adverse pressure 
gradient region. In the region of favorable pressure 
gradient, the TKE model velocities overshoot the PML model 
velocities at some locations and undershoot them at others. 

The above observations can be explained by the fact 
that the mixing length model implies the diffusion and con- 
vection terms of the turbulence kinetic energy equation are 
negligible [14] . One can see this by setting the convection 
and diffusion terms equal to zero in Equation (10) , the 
turbulence kinetic energy equation then becomes : 

u*^ 

L l2i*= 

p 9z* k(z*+ z*) 

where , 


uj = /t/p 
Thus, 

f3u* 
9z* 


T = pJl' 


9u* 

9z* 
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Figure 2. Velocity Distribution Over a 2/1 Ellipse for 
Turbulence Kinetic Energy and. Mixing Length 
Models . 






Figure 3. Velocity Distribution Over a 4/1 Ellipse for 
Turbulence Kinetic Energy and Mixing Length 
Models . 
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Figure 4. Velocity Distribution Over a 10/1 Ellipse 

for Turbulence Kinetic Energy and Mixing 
Length Models. 
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where 


Z*= K ( Z* + Z* ) 

which is the Prandtl mixing length hypothesis. 

The production, diffusion, dissipation, and convec- 
tion terms are plotted versus x in Figures 5, 6, and 7 at 
elevations of z equal to 0.063, 0.206, and 0.42 respectively. 
Also, plotted in each figure are the velocities at the given 
z for the TKE model and the PML model. 

For all three values of z, the production term is 
larger than the dissipation term in the major part of the 
adverse pressure gradient region x < -5. This reflects an 
imbalance between production and dissipation in the adverse 
pressure gradient region with the TKE model which would not 
occur with the PML model, since the production and dissipa- 
tion are equal in the latter model. The production term is 
a measure of the rate at which kinetic energy is extracted 
from the mean flow and converted into turbulence kinetic 
energy. Therefore, more kinetic energy is removed from the 
mean flow in the adverse pressure gradient region for the 
TKE model than for the PML model, resulting in a greater 
deceleration of the flow. 

In the favorable pressure gradient region, x > -5, 
the behavior of the models varies with respect to each other 
depending on the height z. Near the surface (z 0.063), 
Figure 5, the production term of the TKE model is less than 
the dissipation term in the initial portion of the favorable 


32 



Figure 5. Turbulence Kinetic Energy Production 

Dissipation, Diffusion, and. Convection 
Versus Distance Upstream from Top of 
Ellipse at z = 0.063. 
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Figure 6. Turbulence Kinetic Energy Production, Dissipation, 
Diffusion, and Convection Versus Distance Upstream 
from Top of Ellipse at z = 0.206. 
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Figure 7, Turbulence Kinetic Energy Production, Dissipation, 

and Convection Versus Distance Upstream 
from Top of Ellipse at z = 0.42. 
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pressure gradient region (-5 < x < -3 ) , implying that less 
kinetic energy is removed from the mean flow by the TKE 
model than by the PML model. This results in greater 
acceleration of the flow by the pressure gradient so that 
the velocity for the TKE model catches up to and overshoots 
that of the PML model. 

In the portion of the region of favorable pressure 
gradient from x > -3 to x = 0, production again exceeds 
dissipation for the TKE model. More kinetic energy is 
extracted from the mean flow so that the flow predicted by 
the TKE model now accelerates at a lesser rate, allowing the 
velocity for the PML model to overtake the TKE model 
velocity. 

Further from the surface at z = 0.206 and 0.42, 

Figures 6 and 7, production is less than dissipation for the 
TKE model, resulting in greater acceleration by the pressure 
gradient so that the velocities for TKE model overshoot 
those for the PML model at the top of the ellipse. 

The differences in behavior of the production and 
dissipation terms at different heights above the surface can 
possibly be explained by wall effects on the flow very near 
the surface. 

In the region of adverse pressure gradient, the 
combination of friction and pressure gradient near the sur- 
face decelerate the flow there at a greater rate than further 
from the surface where frictional effects are not as large 
and the momentum of the flow attempts to hold the velocity 
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constant. This tends to form a region of higher shear 
analogous to a free shear layer. Experimental evidence of 
this is shown in Figure 8, from [15], which shows peak u'w* 
values between the surface layer and the outer layer of the 
flow. As a result of this layer of increased shear, the 
production term, (t/p) 3u/3z, tends to decrease near the 
surface and increase further away as the stagnation point is 
approached . 

In the favorable pressure gradient region, the 
pressure gradient works against the wall friction which 
causes the layer of maximum shear to move toward the surface 
with the velocity profile becoming flatter further away. 

Thus the production of turbulence kinetic energy increases 
near the surface and decreases at higher elevations. 
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CHAPTER VI- 


EFFECTS OF ASPECT RATIO 

The effects of a variation of elliptical aspect ratio 
on the velocity profiles at the top of the ellipse are 
shown in Figures 9 and 10 for the turbulence kinetic energy 
model and Prandtl mixing length model respectively. For 
both models the effect of increasing the aspect ratio 
results in a decreased overshoot in velocity at the top of 
the ellipse. 

This effect can be explained by the variation with 
aspect ratio of the pressure gradient distribution along the 
x-axis as shown in Figure 11. The regions of adverse 
pressure gradient are very similar for the three bodies. 
Therefore, the retardation of the flow in the adverse pres- 
sure gradient regions will be essentially the same. The 
regions of favorable pressure gradient, however, are quite 
different for the three bodies. The magnitude of the 
favorable pressure gradient increases with decreasing aspect 
ratio. For the 2/1 ellipse, the pressure gradient varies 
markedly over the entire surface of the ellipse whereas for 
the 10/1 ellipse the variation is only over a small region 
of the upstream surface. The flow is therefore much more 
strongly accelerated by the favorable pressure gradient of 
the lower aspect ratio ellipse and as a result, attains 
higher velocities at the top. 
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Figure 9. Velocity Profiles at the Top of the Ellipse for 
Turbulence Kinetic Energy Model for Aspect 
Ratios of 2/1, 4/1, 10/1. 
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Velocity Profiles at the Top of 
Prandtl Mixing Length Model for 
Ratios of 2/1, 4/1, 10/1. 
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Figure 11. Variation of Pressure Gradient Distribution 
Along x-Axis with Respect to Aspect Ratio 
of Ellipse. 
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CHAPTER VII 

INFLUENCE OF SURFACE ROUGHNESS AND ASPECT 
RATIO ON SEPARATION 

The adverse pressure gradients upstream of the body 
cause separation of the boundary layer near the stagnation 
region. The point at which this separation occurs can be 
determined by solving the governing equations discussed 
previously, using the pressure distribution given by the 
potential flow solution for the streamline i)> = 0. This was 
done for a number of surface roughnesses and aspect ratios. 
The results of the calculations are shown in Figure 12 . 

The length of the forward separation bubble, X, 
decreases with decreasing surface roughness and/or increasing 
aspect ratio. This behavior of the separation bubble is 
similar to that for the Prandtl mixing length model [4] , 
except that separation occurs further upstream for the TKE 
model in all cases. This is to be expected because of the 
greater retardation of the flow in the adverse pressure 
gradient region described previously for the TKE model. 
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Effects of Surface Roughness and Aspect Ratio on Length 
of Separation Bubble, A.. 


Figure 12. 



CHAPTER VIII 


TURBULENCE INTENSITY DISTRIBUTION NEAR 
AN ELLIPTICAL OBSTRUCTION 

If an aircraft is to fly over an obstruction in the 
wind field at low altitude, it is necessary to be able to 
predict regions of high turbulence levels so that the pilot 
can avoid these areas or at least be prepared to encounter 
the turbulence. 

The distribution of turbulence intensity, 

fi' = *^1/3 (u"' ^ + V ' 2 + w ' 2 ) /U^ (34) 

over an elliptical body with an aspect ratio of 4/1 is shown 
in Figure 13. The turbulence intensity is greatest near the 
body and is nearly uniform above an altitude of twice the 
height of the ellipse. There is a local maximum in turbu- 
lence intensity above the stagnation point. 

The turbulence intensity distribution is proportional 
to the square root of the turbulence kinetic energy distri- 
bution, since, from Equations (2) and (34) 

a* = /2/3 e/U^ (35) 

The turbulence kinetic energy distribution is shown in 
Figure 14. The variations of the turbulence kinetic energy 
in the flow field is related to the behavior of the convec- 
tion term of the turbulence kinetic energy equation (Figure 
15) . The convection term indicates the rate of convection 
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of turbulence kinetic energy toward or away from a given 
point. A positive value of the convection term in a given 
area indicates convection away from the area, while a 
negative value indicates convection toward the area. The 
areas where the turbulence kinetic energy is increasing in 
the downstream direction correspond to areas where the con- 
vection term is positive, and areas of decreasing turbulence 
kinetic energy correspond to a negative convection term. 

The convection term is the sum of the production and 
diffusion terms (Figures 16 and 17 respectively) minus the 
dissipation term (Figure 18) . The sign of the convection 
term is primarily determined by the production and dissipa- 
tion terms, which are an order of magnitude greater than the 
diffusion term in most areas. A positive convection term 
indicates that production is greater than dissipation, 
requiring convection of turbulence kinetic energy to sur- 
rounding areas, the reverse being true in areas of negative 
convection. 

The peak in turbulence kinetic energy near the stag- 
nation point occurs where the flow near the surface is 
retarded by the friction at the wall, while the flow at a 
greater distance from the wall is accelerated by a favorable 
pressure gradient, which decreases with increasing height to 
match the free stream logarithmic velocity profile. The 
steep turbulence kinetic energy (and turbulence intensity) 
gradients occurring near the wall in this area are due to 
the transition from an adverse to a favorable pressure 
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Figure 17. Distribution of Turbulence 
4/1 Ellipse. 



Kinetic Energy Diffusion Over 
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gradient . 

The slight local minimum in turbulence kinetic energy 
between x -1 and x = 0 at z 2 is due to the velocity 
decreasing from its peak value at z ^ 1.5 in order to match 
the free stream logarithmic profile (Figure 3/ page 30). 

This results in a region of low shear stress near z = 2. 

An aircraft flying over this elliptical body at an 
altitude less than twice the ellipse height would experience 
increased turbulence with decreased wind velocity at 
X -4, and increased turbulence and wind velocity over the 
top of the ellipse. The intense turbulence zones may prove 
to be unexpected hazards during landing over elliptical 
surface obstructions. 
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CHAPTER IX 


DISCUSSION OF RESULTS RELATIVE TO 
AERONAUTICAL APPLICATIONS 

An assessment of the influence of the wind disturbance 
created by a semi-elliptical terrain irregularity on air- 
craft operations through the disturbance is given in this 
section. The wind environment encountered along the flight 
paths for CTOL and STOL aircraft relative to the FAA surface 
obstruction clearance recommendations is discussed in rela- 
tion to experiments and analytical studies reported in the 
literature on aircraft airworthiness as influenced by wind. 

Turbulence 

If one considers the elliptical surface to be a hill 
over which approaches are made, the flight paths for CTOL 
and STOL aircraft based on FAA surface obstruction clearance 
standards (see [16, 17]) would be as shown in Figure 19. In 
general the aircraft passes through only the fringes of the 
disturbance with the STOL' aircraft experiencing at most a 
region of turbulence kinetic energy approximately 10 percent 
higher than the undisturbed atmospheric value. Assuming a 
9 m/sec wind in the undisturbed flow at a height of three 
times that of the hill, gives = 1.41 m/s in the undis- 
turbed region and = 1.48 m/s in the portion of the dis- 
turbed region through which the airplane passes. 

Neuman and Foster [18] have observed that with 
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= 1.83 the tentative FAA spec, fication of 0.9 m/s maximum 
rate of descent is exceeded 45 pe. cent of the time. 

Assuming two Gaussian distributions with standard deviations 
of 1.41 m/s and 1.48 m/s, respectively, the probability of 
exceeding 1.83 m/s for each distribution can be predicted. 

The probability of a hard landing is then 45 percent , times 
the probabili-ty that 1.83 m/s will be exceeded. Based on 
this analysis there is a 10 percent higher chance of a hard 
landing than in the natural atmosphere which due to the crude 
nature of the analysis is seemingly negligible. Hence the 
FAA clearance surface requirements for STOL ports appear 
adequate for approaches over the windward part of curved 
proturbances . However, if one envisions landing on a strip 
on top of a plateau or hill, the aircraft will encounter 
considerably higher turbulence. On the top of the elliptical 
surface, the ratio of the rms components in the disturbed 
flow to that in the undisturbed flow can be as high as 1.58 
given = 2.23 m/s. The probability of hard landings in 
this case is, according to the crude analysis, 100 percent 
greater than in the undisturbed atmosphere. The above con- 
clusions are obviously overly pessimistic since the more 
intense turbulence considered here occurs only in patches 
whereas the results of Nexaman and Foster [18] assume con- 
tinuous homogenous turbulence. 

Mean Wind 

Figure 20 shows the variation in the horizontal mean 
wind along the flight path for STOL and CTOL aircraft over 
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the 4/1 and 2/1 ellipse respectively. For the STOL aircraft 
passing over the 2/1 ellipse there is approximately a 2/ms 
change in velocity over a horizontal distance of eight semi- 
ellipse heights. The CTOL aircraft "sees" less wind varia- 
tion along the glide slope due to its less steep approach 
path. If one considers an airplance landing at 80 kts over 
a 5 m hump the 2 m/s (approximately 4 kts) velocity change 
will be experienced in approximately 0.97 sec. Adding the 
= 1.48 m/s value to the 2 m/s value the variation in wind 
will appear to the aircraft as a gust approaching the order 
of the 6 m/s lateral gust value over a 0.7 sec period 
reported by McManus [1] . These gusts were observed during 
experiments with automatic landing of the Trident and their 
effect was that the autopilot controlled the bank during the 
increase in gust but could not keep pace with the decay, 
consequently the aileron drove the aircraft away from the 
wings level condition. Although not disastrous, such gusts 
did cause the aircraft to deviate significantly from the 
desired flight path while the autopilot reestablished 
integration. 

In this same vein, Corbin and Goddard [19] have 
examined the effect of longitudinal gusts on experimental, 
automatic flight control systems for the BAG 1-11. Their 
predicted maximum rate of descent at touchdown as a function 
of gust strength is shown in Figure 21. Gusts of 2 m/s (4 kt) 
magnitude are observed to give touchdown rates between 1.1 
to 0.6 m/s which brackets the tentative FAA maximum descent 
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rate of 0.9 m/s [20]. 

Wind Shear 

The local wind shear encountered along the STOL 
flight path is shown in Figure 22. Severe gradients which 
change sign near the front of the ellipse are observed. For 
comparison purposes, the shear along a 6 degree glide slope 
in a wind over natural, homogenous terrain (i.e., in the up- 
stream, logarithmic wind profile) is also shown. 

The influence of spatially varying wind shears on 
aircraft stability has not been investigated in general 
form, however, studies of uniform wind shear (linear wind 
profiles) are reported in [21, 22, 2]. Figure 23, from 
Etkin [21] shows the effect of constant wind shear on 
phugoid roots and short-period roots for a STOL airplane. 

The range of wind shear from +0.30 to -0.30 sec ^ is covered 
in the figure, which it should be noted, does not encompass 
the full +0.40 to -0.50 sec ^ range which would be 
encountered passing over the semi-ellipse. 

The effect of wind shear on the phugoid and pitching 
modes are seen to be large. A strong positive shear 
decreases both the frequency and damping of the phugoid, and 
a strong negative shear changes the real pair of pitching 
roots into a complex pair representing a pitching oscilla- 
tion of long period and heavy damping. The time the air- 
craft is in the disturbance is short compared to the phugoid 
period and hence the phugoid oscillation is unimportant. 
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Figure 22. Wind Shear Encountered By STOL Aircraft 




However, the pitching oscillations would be significant. 
These effects could become large in the disturbed wind over 
the semi-ellipse. It should be noted that the analyses of 
[21] assumes the aircraft continuously experiences shear, 
whereas the shear discussed here is highly localized. 

Gera [22] carried out a similar analysis to Etkin [21] 
and found that for his particular airplane an unstable root 
appears for wind shears of 0.183 sec"^. He concludes that 
for certain configurations and flight conditions the appear- 
ance of this diverging first-order root may be important and 
may warrant careful analysis. 

Conclusions 

The preceding section has described the computed flow 
fields in context with aircraft operations. Comparisons are 
made with other reported analyses and experiments based on 
idealized models of both the flow fields and the airplane 
dynamics, however, the following conclusions are warranted. 

In all the comparisons made the magnitude of the winds com- 
puted in this study along the 6 degree STOL flight path are 
on the order of wind conditions considered adverse by other 
investigators. This is not the case for CTOL aircraft 
because of the less steep flight path, 3 degrees. 

The recommended 15:1 obstruction clearance surfaces 
[17] for STOL ports in urban areas, therefore, appear per- 
haps questionable in view of the fact that wind disturbances 
created by semi-elliptical terrain regularities extend into 
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the clear region sufficiently far to generate wind environ- 
ments which are boarder-line for safe operations. Addi- 
tionally bluff bodies and trailing wakes behind surface 
proturbances will create even stronger wind disturbance than 
the smooth curved ellipse, although, three dimensionality 
will tend in some cases to alleviate these disturbances. 

No attempt has been made here to conduct an extensive 
parametric study of wind speeds, surface roughness, aspect 
ratio, etc. on the airplane performance. Rather an effec- 
tive tool embodied in the computer code for modelling the 
disturbed wind field is presented which can be used to 
conduct a more generalized study. In view of the apparent 
closeness of the wind environment over an elliptical surface 
to boarder-line hazardous flight winds, careful appraisal of 
current thinking on 6 degree glide slopes and 15:1 obstruc- 
tion clearance surfaces seems necessary and further investi- 
gation along these lines is being continued. 
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CHAPTER X 


CONCLUSIONS 

The boundary layer analyses presented herein of flow 
over a semi-elliptical surface obstruction using the turbu- 
lence kinetic energy equation leads to the following 
conclusions : 

1. The Prandtl mixing length model and turbulence 
kinetic energy equation model predict essentially 
the same velocity profiles for the boundary layer 
analysis of flow over a semi-elliptical surface. 
Slight variations in the wind profiles occur in 
the region of varying pressure gradient due to 
the diffusion and convection properties of the 
turbulence kinetic energy inherent to the TKE 
model and absent in the PML model. It is believed 
for this reason the TKE model gives somewhat 
better velocity prediction than the PML model. 

2 . A region of high turbulence intensity occurs 
above the stagnation point at approximately the 
height of the semi-ellipse and a sharp gradient 

in intensity appears along a line where transition 
from adverse to favorable pressure gradient takes 
place . 

3. The same conclusions are made regarding the 
influence of surface roughness and ellipse aspect 
ratio on maxima in wind speed upstream separation 
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bubble and velocity profiles as were drawn for 
the PML model and stated earlier on page 3 of 
this paper. 

4. The code developed in this study provides an 

effective tool for investigating the influences 
of a wind disturbed by a terrain irregularity on 
the flight operations of aircraft. Interpreta- 
tion of the computed flow fields over the semi- 
elliptical surface suggests that the 15:1 
obstruction clearance and 6 degree glide slope 
presently recommended for STOL aircraft may be 
inadequate for saft flight operations. 
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